%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This program replicates Figure 3 in the paper "Optimal 
%%% Portfolio Choice with Fat Tails and Parameter Uncertainty",
%%% Journal of Financial and Quantitative Analysis (2024),
%%% Raymond Kan & Nathan Lassance 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
close all

%%% Exponential integral
Expfun = @(n,x) integral(@(t)exp(-x.*t)./(t.^n),1,inf);

%%% Initialize vector
rhovec = 0:0.01:1;
ratioc1psi02nu4 = zeros(length(rhovec),1);
ratioc1psi04nu4 = zeros(length(rhovec),1);
ratioc2psi02nu4 = zeros(length(rhovec),1);
ratioc2psi04nu4 = zeros(length(rhovec),1);
ratioc1psi02nu8 = zeros(length(rhovec),1);
ratioc1psi04nu8 = zeros(length(rhovec),1);
ratioc2psi02nu8 = zeros(length(rhovec),1);
ratioc2psi04nu8 = zeros(length(rhovec),1);

%%% \nu = 4
nu = 4;
j = 1;
for rho = rhovec
    if rho == 0
        ratioc1psi02nu4(j,1) = 1;
        ratioc1psi04nu4(j,1) = 1;
        ratioc2psi02nu4(j,1) = 1;
        ratioc2psi04nu4(j,1) = 1;
    elseif rho == 1
        %\psi = 0.2
        psi2 = 0.2^2;
        ratioc1psi02nu4(j,1) = (psi2+rho)/((nu/(nu-2))*psi2+rho);
        ratioc2psi02nu4(j,1) = ((nu-2)/nu)*(psi2+rho)/((nu/(nu-2))*psi2+rho);
        %\psi = 0.4
        psi2 = 0.4^2;
        ratioc1psi04nu4(j,1) = (psi2+rho)/((nu/(nu-2))*psi2+rho);
        ratioc2psi04nu4(j,1) = ((nu-2)/nu)*(psi2+rho)/((nu/(nu-2))*psi2+rho);
    else
        y = @(x) (nu-2).*rho.*x./(2*(1-rho));
        fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
        eta = fzero(fun,0.5*(1+nu/(nu-2)));
        varphi = 2*eta^2*(1-rho)/(nu-eta*(nu-2));
        %\psi = 0.2
        psi2 = 0.2^2;
        ratioc1psi02nu4(j,1) = (psi2+rho)/(psi2*varphi/eta+rho);
        ratioc2psi02nu4(j,1) = (eta/varphi)*(psi2+rho)/(psi2*varphi/eta+rho);
        %\psi = 0.4
        psi2 = 0.4^2;
        ratioc1psi04nu4(j,1) = (psi2+rho)/(psi2*varphi/eta+rho);
        ratioc2psi04nu4(j,1) = (eta/varphi)*(psi2+rho)/(psi2*varphi/eta+rho);
    end
    j = j+1;
end

%%% \nu = 8
nu = 8;
j = 1;
for rho = rhovec
    if rho == 0
        ratioc1psi02nu8(j,1) = 1;
        ratioc1psi04nu8(j,1) = 1;
        ratioc2psi02nu8(j,1) = 1;
        ratioc2psi04nu8(j,1) = 1;
    elseif rho == 1
        %\psi = 0.2
        psi2 = 0.2^2;
        ratioc1psi02nu8(j,1) = (psi2+rho)/((nu/(nu-2))*psi2+rho);
        ratioc2psi02nu8(j,1) = ((nu-2)/nu)*(psi2+rho)/((nu/(nu-2))*psi2+rho);
        %\psi = 0.4
        psi2 = 0.4^2;
        ratioc1psi04nu8(j,1) = (psi2+rho)/((nu/(nu-2))*psi2+rho);
        ratioc2psi04nu8(j,1) = ((nu-2)/nu)*(psi2+rho)/((nu/(nu-2))*psi2+rho);
    else
        y = @(x) (nu-2).*rho.*x./(2*(1-rho));
        fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
        eta = fzero(fun,0.5*(1+nu/(nu-2)));
        varphi = 2*eta^2*(1-rho)/(nu-eta*(nu-2));
        %\psi = 0.2
        psi2 = 0.2^2;
        ratioc1psi02nu8(j,1) = (psi2+rho)/(psi2*varphi/eta+rho);
        ratioc2psi02nu8(j,1) = (eta/varphi)*(psi2+rho)/(psi2*varphi/eta+rho);
        %\psi = 0.4
        psi2 = 0.4^2;
        ratioc1psi04nu8(j,1) = (psi2+rho)/(psi2*varphi/eta+rho);
        ratioc2psi04nu8(j,1) = (eta/varphi)*(psi2+rho)/(psi2*varphi/eta+rho);
    end
    j = j+1;
end

%%% Create Figure 3
fig = figure;
%Subplot 1
subplot(2,2,1);
plot(rhovec,ratioc1psi02nu4,'b','LineWidth',2);
hold on
plot(rhovec,ratioc1psi04nu4,'-.b','LineWidth',2);
plot(rhovec,ratioc1psi02nu8,'--r','LineWidth',2);
plot(rhovec,ratioc1psi04nu8,':r','LineWidth',2);
hold off
legend('$(\nu,\psi)=(4,0.2)$','$(\nu,\psi)=(4,0.4)$','$(\nu,\psi)=(8,0.2)$',...
       '$(\nu,\psi)=(4,0.4)$','Location','southeast','interpreter','latex');
set(gca,'TickLabelInterpreter','latex');
grid on
title('$c_1^\star/c_{1,n}^\star$','interpreter','latex');
xlabel('$\rho$','interpreter','latex');
xlim([0 1]);
ylim([0.8 1]);
xticks([0 0.2 0.4 0.6 0.8 1]);
yticks([0.8 0.85 0.9 0.95 1]);
%Subplot 2
subplot(2,2,2);
plot(rhovec,ratioc2psi02nu4,'b','LineWidth',2);
hold on
plot(rhovec,ratioc2psi04nu4,'-.b','LineWidth',2);
plot(rhovec,ratioc2psi02nu8,'--r','LineWidth',2);
plot(rhovec,ratioc2psi04nu8,':r','LineWidth',2);
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('$c_2^\star/c_{2,n}^\star$','interpreter','latex');
xlabel('$\rho$','interpreter','latex');
xlim([0 1]);
ylim([0.4 1]);
xticks([0 0.2 0.4 0.6 0.8 1]);
yticks([0.4 0.6 0.8 1]);
fontsize(fig,16,"points");